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I.  INTRODUCTION 
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Proton  transfer  in  water  has  been  the  focus  of  considerable  theoretical  and  experimental 
investigation  during  this  century.  A  lone  proton  does  not  exist  in  water  but  instead  forms  a 
weak  chemical  bond  with  a  water  molecule  to  give  a  hydronium  (  HsO'*'  )  ion.  The  mobility 
of  the  hydronium  itself  is  anomalously  high,  being  five  to  seven  times  that  of  similarly  sized 
cations  [1].  The  explanation  of  this  unusual  fact  is  thought  to  be  that  the  migration  of  an 
excess  proton  in  water  could  occur  via  a  mechanism  of  chemical,  rather  than  hydrodynamic 
or  Stokes’  law,  diffusion.  The  positive  charge  might  therefore  shuttle  through  water  via 
an  exchange  of  chemical  and  hydrogen  bonds  (as  illustrated  in  Fig.  1),  hence  the  charge 
migration  could  occur  via  a  process  of  structural  diffusion,  [2]  also  called  the  Grotthus 
mechanism. 

A  classic  variant  on  this  model  suggested  by  Bockris  [3,4]  is  one  in  which  the  rate  deter¬ 
mining  step  is  the  molecular  reorientation  of  a  water  adjacent  to  the  in  order  to  form 

a  hydrogen  bond  with  the  HsO"*”  ion.  One  drawback  of  this  theory  is  that  the  thermal  rate 
of  reorientation  of  water  molecules  in  pure  water  is  insufficiently  fast  to  account  for  the  rate 
of  proton  transfer  [3].  However,  it  can  be  postulated  that  the  electric  field  of  the  hydronium 
enhances  the  rate  of  molecular  reorientation.  This  process  is  known  as  “field  assisted”  re¬ 
orientation  and  is  could  be  sufficiently  fast  to  account  for  the  rate  of  proton  transfer  [4] .  On 
the  other  hand,  it  has  been  pointed  out  [5-7]  that  the  activation  energy  necessary  to  transfer 
a  proton  to  an  acceptor  water  that  is  not  hydrogen  bonded  to  other  water  molecules  in  the 
second  solvation  shell  (which  is  a  consequence  of  the  field  reorientation  mechanism)  exceeds 
the  experimentally  determined  [8]  activation  energy  of  proton  transfer.  Furthermore,  the 
theory  as  originally  put  forward  [3,4]  would  seem  to  require  the  presence  of  a  second  hydro¬ 
nium  to  perturb  the  local  structure  around  the  first  hydronium  (see  also  [9]).  In  addition, 
the  mobilities  of  H+  and  OH“  are  smaller  [10]  in  ice  than  in  water  which  would  seem  to 
contradict  a  field  “assisted”  reorientation  mechanism.  Hence,  it  seems  that  the  rate  limit¬ 
ing  step  for  proton  transfer  might  not  involve  reorientation  dynamics  in  the  first  solvation 
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shell  but  instead  could  involve  dynamics  of  water  molecules  in  the  second  solvation  shell  of 
the  HaO"*"  ion.  In  this  picture,  the  formation  of  hydrogen  bonds  in  the  second  solvation 
shell  with  a  water  molecule  in  the  first  solvation  shell  is  a  precursor  to  a  Grotthus— type 
proton  transfer  event  [7,11-13].  Also  critical  to  a  successful  proton  transfer  event  could  be 
the  inward  fluctuation  of  the  oxygen-oxygen  distance  between  the  hydronium  ion  and  the 
acceptor  water  molecule  in  the  first  solvation  shell  [14,15]  since  the  barrier  to  proton  transfer 
increases  rapidly  as  the  oxygen-oxygen  distance  increases  [16-20].  Hence,  proton  transfer 
must  be  a  concerted  process  depending  on  fluctuations  in  both  the  oxygen-oxygen  distance 
between  donor  and  acceptor  as  well  as  the  nearby  solvent  dynamics. 

The  proton  due  to  its  small  mass  can  give  rise  to  quantum  effects  in  many  chemical 
reactions  (see,  e.g..  Refs.  [21-25]).  The  thermal  DeBroglie  wavelength  of  a  proton  is  ~  1.0 
A  which  is  of  the  same  magnitude  as  the  length  scale  over  which  proton  transfer  takes  place 
in  solution.  In  addition,  the  isotope  ratios  for  the  rate  constants  of  the  chemical  reactions 
H2O  +  H3O+  H3O+  +  H2O  and  H2O  +  OH'  OH"  +  H2O  have  been  measured 
[26]  to  be  kf  /  kf  =  1.6  ±  0.2  and  kf  /  k^  =  2.7  ±  0.4,  respectively.  The  effects  are  somewhat 
larger  than  the  kinetic  isotope  effect  which  indicates  that  quantum  effects  may 

be  present  to  some  degree  in  both  reactions. 

The  hydronium  species  [4,11,12]  in  water  has,  on  the  average,  €3^  symmetry  [27]  and 
four  water  molecules  in  its  first  solvation  shell  [28].  Alternatively,  Zundel  [29]  and  others  [30] 
have  characterized  the  excess  proton  as  an  HsO^  “grouping”  in  water  in  which  the  excess 
proton  tunnels  rapidly  back  and  forth  between  the  two  water  molecules  through  the  strong 
hydrogen  bond.  There  is  some  contention  in  the  literature  as  to  which  species  predominates 
in  the  liquid  [27].  A  recent  pioneering  ab  initio  MD  study  [31]  indicates  that  both  viewpoints 
are  in  essence  correct,  at  least  for  classical  nuclear  dynamics.  In  this  study,  it  was  found 
that  there  is  a  rapid  interconversion  between  the  H9O4  and  HsO^  structures  [32]  in  solution 
and  this  interconversion  is  governed  by  the  dynamics  of  the  local  solvent  structure  in  the 
second  solvation  shell  of  the  H30'*'  ion. 

Ab  initio  studies  of  proton  transfer  in  small  gas  phase  water  clusters  are  numerous  (see. 
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e.g.,  Refs.  [6,7,16-20,33-37]).  These  studies  have  characterized  the  H5OJ  complex  as  having 
a  flat  and  broad  potential  on  which  the  excess  proton  moves.  The  proton  potential  surface, 
depending  on  the  level  of  ab  initio  theory  used,  is  either  a  single  minimum  or  double  minimum 
with  a  transition  state  0.2-0.5  kcal  mol"^  higher  than  the  minima.  As  additional  water 
molecules  are  added,  the  potential  of  the  proton  becomes  a  single  minimum  and  shifts 
according  to  the  number  and  positioning  of  water  molecules  added  to  the  complex  [6,7].  In 
extending  this  picture  to  the  bulk  solvent,  it  is  thought  the  dynamics  of  the  excess  proton 
is  governed  by  solvent  fluctuations  about  the  proton  transfer  complex.  The  electric  field 
of  the  solvent  determines,  in  part,  the  potential  on  which  the  proton  moves  and  induces 
shifts  of  the  proton  from  one  water  molecule  to  the  other.  As  the  proton  shifts  there  is  a 
charge  redistribution  in  the  proton  transfer  complex  (  H5OJ  ).  This  charge  redistribution  is 
a  polarization  effect  termed  “Zundel  Polarization”  [38-40]. 

There  have  been  few  computational  investigations  of  proton  transfer  in  water  using 
Molecular  Dynamics  (MD)  or  Monte  Carlo  (MC)  methods  [41-44]  because  of  the  difiiculty 
in  specifying  empirical  potential  functions  that  correctly  describe  proton  transfer  along  a 
series  of  hydrogen  bonded  water  molecules.  The  recent  “first-principles”  studies  of  this 
and  related  systems  [31,45-48]  using  ab  initio  MD  [also  known  as  the  Car-Parrinello  (CP) 
method  [49-52]]  have  taken  a  significant  step  forward  forward  on  the  problem.  The  CP 
method  is  a  dynamical  density  functional  theory-based  (DPT)  method  wherein  the  elec¬ 
tronic  structure  is  calculated  “on  the  fly”  while  being  adiabatically  “slaved”  to  the  evolving 
classical  nuclei  within  the  Born-Oppenheimer  approximation.  The  forces  on  the  classical 
nuclei  are  calculated  from  the  Hellman-Feynman  theorem.  Hence,  chemical  processes  such 
as  bond  breaking  or  formation  are  treated  at  an  ab  initio  level  rather  than  using  empirical 
methods  and/or  potentials,  representing  a  substantial  advance  in  MD  simulation  capabil¬ 
ity.  Unfortunately,  the  high  computational  cost  has  limited  the  size  of  systems  studied  and 
prevented  the  length  of  the  CP  simulations  from  being  more  than  four  to  five  picoseconds, 
which  is  often  insufficient  to  obtain  well  converged  values  of  thermodynamic  and  dynami¬ 
cal  quantities  of  interest.  Moreover,  the  accuracy  of  gradient-corrected  DPT  methods  for 


4 


describing  the  proton  transfer  barrier  in  the  H5O2  complex  has  been  called  into  question 
[15],  and  the  nuclei  are  treated  as  being  classical  in  the  CP  simulations. 

In  the  present  paper,  the  dynamical  and  equilibrium  properties  of  an  excess  proton 
are  studied  via  classical  MD  and  Feynman  path  integral  methods.  The  potential  energy 
modeling  is  accomplished  using  a  two  state  Empirical  Valence  Bond  (EVB)  [53,54]  model 
for  the  H5OJ  complex  in  water.  The  goal  of  this  work  is  somewhat  more  modest  than 
an  explicit  simulation  of  the  Grotthus  mechanism:  the  activation  free  energy  as  well  as  the 
dynamics  of  a  proton  transferring  between  a  hydronium  ion  and  an  acceptor  water  molecule 
in  the  first  solvation  shell  are  calculated  using  both  the  quantum  and  classical  treatments  of 
the  system.  The  quantum  effects  on  the  proton  exchange  along  the  strong  hydrogen  bond 
between  the  H3O+  ion  and  a  water  molecule  are  thus  examined.  From  the  results  presented 
in  the  following  sections,  it  will  be  suggested  that  the  transfer  of  a  proton  from  a  hydronium 
ion  to  a  neighboring  water  molecule  is  not  the  rate  limiting  step  in  the  proton  migration 
process.  Rather,  the  likely  quantum  dynamical  mechanism  for  extended  (Gr5tthus-type) 
proton  transfer  can  be  uncovered  by  examining  the  structure  and  dynamics  of  the  solvent 
molecules  around  the  proton  transfer  complex.  To  our  knowledge,  the  present  work  is  the 
first  quantum  simulation  of  the  nuclear  dynamics  of  an  excess  proton  in  water. 

The  sections  of  this  paper  are  organized  as  follows:  In  Sec.  II,  the  computational  methods 
used  to  simulate  the  dynamics  and  energetics  of  an  excess  proton  in  water  are  outlined.  The 
parameterization  of  the  empirical  model  is  then  described  in  Sec.  Ill,  while  the  results  of  the 
quantum  and  classical  simulations  are  discussed  in  Sec.  IV .  Concluding  remarks  are  given 
in  Sec.  V. 
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II.  COMPUTATIONAL  METHODS 
A.  Path  Integral  Quantum  Transition  State  Theory 

In  several  previous  studies  of  proton  transfer  reactions  [55-58]  the  Feynman  path  integral 
quantum  transition  state  theory  (PI-QTST)  approach  [59-62]  has  been  used  to  calculate 
the  quantum  mechanical  activation  free  energy  of  proton  transfer.  Both  solvent  and  in¬ 
tramolecular  contributions  to  the  activation  free  energy  of  proton  transfer  are  included  in 
this  approach  and  the  method  is  applicable  [63]  in  both  the  adiabatic  and  nonadiabatic 
limits  of  proton  transfer  and  where  the  proton  itself  becomes  classically  activated.  Impor¬ 
tantly,  the  quantum  tunneling  of  the  proton  is  explicitly  included  in  the  calculation.  [56,63] 
The  PI-QTST  method  makes  no  approximation  or  assumptions  concerning  the  form  of  the 
couphng  of  the  proton  transfer  coordinate  to  the  solvent  and  includes  the  highly  nonlinear 
intramolecular  and  intermolecular  couplings  between  the  proton  transfer  coordinate  and  the 
other  molecular  modes  of  the  proton  transfer  complex. 

In  the  PI-QTST  formulation,  [59-61]  the  fundamental  quantum  rate  constant  for  a 
quantum  activated  rate  process  can  be  written  as 

^  ,  (2.1) 

hQR 

where  k  is  a  quantum  dynamical  transmission  coefficient  [59-61]  of  order  unity,  Qr  is  the 
partition  function  of  the  entire  solute-solvent  system  in  its  reactant  configuration,  ks  is 
Boltzmann’s  constant,  T  is  the  temperature,  /3  is  l/ksT,  and  F*  is  a  quantum  mechanical 
excess  free  energy,  given  by  [59-61] 

F;  =  -kBTln[Q,{q*)/{f^/27rn^0Y^^]  .  (2-2) 

The  quantity  Qdq*)  in  Eq.  2.2  is  the  path  integral  centroid  density  [59-61,64],  given  in  the 
case  of  a  proton  transfer  reaction  by  [56,63] 

Qc{q*)  =  J■■■J'DriT)VR{r)5iq*-qo)exp{-S[v{T),R{r)]/n}  ,  (2.3) 
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where  the  coordinates  r  are  the  Cartesian  coordinates  of  the  proton  and  the  coordinates  R 
are  all  other  general  coordinates  of  the  system,  including  the  solvent.  The  coordinate  q  is 
the  proton  transfer  reaction  coordinate  having  a  reduced  mass  fx.  The  reaction  coordinate  in 
the  case  of  H2O-H-OHJ  complex  is  the  asymmetric  stretch  mode  involving  the  transferring 
proton  and  the  two  oxygen  atoms  of  the  donor  and  acceptor  water  molecules,  i.e., 

qasym  =  -  ^'2)  ,  (2-4) 

where  ri  and  r2  are  the  distances  between  the  proton  (labeled  H3  in  Fig.  2)  and  each  oxygen 
atom  respectively  (labeled  Oi  and  O2  in  Fig.  2).  The  centroid  variable  of  the  reaction 
coordinate  path  qij)  is  defined  as  [59-62,64] 

and  is  constrained  to  be  at  the  transition  state  between  the  reactant  and  product  states  of 
the  proton  [56,63]  (which  in  this  case  is  qa5ym“0).  The  action  5[r(r),R(T)]  in  Eq.  (2.3)  is 
the  imaginary  time  path  integral  action  functional  [64],  which  depends  both  on  the  paths 
of  the  proton  coordinates  r(r)  and  those  of  the  remaining  coordinates  of  the  system  R(t). 
The  specific  expression  for  the  action  functional  in  Eq.  (2.3)  is  given  by  [64] 

5[r(r),R(r)l  =  +  +  V(rW.RWl}  .  (2.6) 

where  y[r(r),R(r)]  is  the  general  many-body  potential.  All  possible  quantum  paths  of 
the  coordinates  are  integrated  over  in  the  functional  integration  in  Eq.  (2.3).  The  classical 
approximation  for  any  of  the  variables  R  in  Eq.  (2.3)  is  to  replace  the  quantum  paths  R(t) 
by  their  equivalent  classical  coordinates.  Due  to  the  heavy  mass  of  the  oxygen  atom,  the 
oxygen  coordinates  were  treated  classically  and  so  the  aforementioned  replacement  was  made 
for  their  quantum  paths.  All  of  the  protons  (not  just  the  transferring  proton)  in  the  solvent 
water  molecules  and  transfer  complex  were  treated  quantum  mechanically  in  this  study. 
Equation  (2.1)  can  be  written  in  a  form  particularly  well  suited  for  computer  simulation 


(2.7) 


k  =  K^exp[-Pi^F') 

ZTT 

where  is  the  curvature  of  the  reaction  coordinate  centroid  free  energy  around  the  reactant 
configurations  (i.c.,  for  qq  near  q^),  the  prefactor  k  is  approximated  to  be  unity,  and  the 
difference  in  centroid  free  energies  between  the  reactant  and  transition  state  is  given  by 


AF:  =  -  kBTln[Qciq*)/Qc{qr)] 

=  -  kBTln[Pc{qr-^q*)]  •  (2-8) 


A  number  of  simulation  methods  are  suitable  for  evaluating  path  integrals  and  Eq.  (2.8). 
In  the  case  of  high  barriers,  one  method  is  to  sample  the  classical  coordinates  and  the  path 
integral  quasiparticle  coordinates  using  Path  Integral  Monte  Carlo  (PIMC)  [65-69]  or  Path 
Integral  Molecular  Dynamics  (PIMD)  [70]  with  umbrella  sampling  [71].  However,  since  the 
activation  energy  of  proton  transfer  is  quite  low  (on  the  order  of  kBT)  in  this  particular 
study,  the  probability  P^qr  q*)  is  readily  calculated  by  binning  the  trajectory  of  the 
centroid  of  the  reaction  coordinate  as  it  evolves  unconstrained  (i.e.  without  the  influence  of 
an  umbrella  potential)  during  the  simulation.  The  PIMD  method  was  chosen  to  simulate 
the  proton  transfer  system  in  this  study.  In  this  method,  the  primitive  representation  of  the 
discretized  path  integral  action  functional  is  used  so  that  each  quantum  particle  maps  onto 
an  isomorphic  polymer  “necklace”  consisting  of  P  quasiparticles,  held  together  by  harmonic 
springs  [72,73].  The  partition  function  in  a  notation  specific  to  the  proton  transfer  complex 


and  the  accompanying  system  of  water  molecules  is  given  by 

I  no0±2  I 

i=l  \j=l 


Qp  — 


/  \  3A/'o/2  p  / \ 

(m)  n  n<  *'‘’n‘'Ro,exp(-/3yp)  (2.9) 


where  the  proton  (oxygen)  coordinates  of  the  water  solvent  as  well  as  donor  and  acceptor 
waters  are  denoted  by  Rh(Ro)*  The  proton  (oxygen)  mass  is  mu  (mo).  The  coordinates  of 
the  transferring  proton  are  denoted  by  r.  The  potential  Vp  in  Eq.  (2.9)  is  given  by 

^  (Rg+b  _  +  Ui+i)  _  rW)' 

+  4  E  ^  ’  (2.10) 


i=l 
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and  consists  of  two  terms.  The  first  term  contains  the  harmonic  springs  for  the  path  inte¬ 
gral  coordinates  which  are  just  the  kinetic  energy  terms  of  Eq.  (2.6).  The  second  term  is 
the  many-body  potential  which  includes  all  the  intermolecular  and  intramolecular  interac¬ 
tions  and  is  evaluated  at  P  discrete  imaginary  time  slices.  The  partition  function  Qp  can 
be  evaluated  as  an  effective  classical  partition  function  by  introducing  fictitious  momenta 
for  the  quasiparticles  and  oxygen  atoms  and  using  molecular  dynamics  [70].  The  effective 
Hamiltonian  is  thus  given  by 


P  /Nn 

EE 

\  3 


(2.11) 


where  is  the  fictitious  mass  of  the  proton  quasiparticles  and  mo  is  the  mass  of  the  oxygen 
nuclei.  The  velocities  of  the  proton  quasiparticles  and  classical  oxygen  nuclei  are  given  by 


and  Vo,  respectively,  while  that  for  the  transferring  proton  quasiparticles  are  given  by 


B.  Proton  Quantum  Dynamics 


The  calculation  of  quantum  dynamical  properties  for  many-body  systems  with  many 
quantum  degrees  of  freedom  is  problematical.  A  number  of  techniques  for  computing  quan¬ 
tum  dynamics  are  available  such  as  semiclassical  methods  [74-80]  or  nonadiabatic  dynamics 
[81-87].  However,  none  of  these  methods  are  particularly  useful  where  a  large  fraction  of  the 
total  nuclei  must  be  quantized.  A  recently  developed  method  known  as  “Centroid  Molecular 
Dynamics”  (CMD)  [88-91]  which  is  capable  of  dealing  with  such  problems  will  be  used  in 
the  present  study  to  simulate  the  dynamics  of  the  proton  transfer  in  a  quantized  water  sol¬ 
vent.  The  focus  of  the  latter  method  is  on  the  path  centroid  as  a  quasi-classical  dynamical 
variable.  The  centroid  of  any  proton  is  defined  in  Eq.  (2.5)  and  in  the  discrete  representation 
is  written  as 


9 


where  Rh(t)  is  the  Cartesian  imaginary  time  paths  of  the  protons  and  are  the  Cartesian 
coordinates  of  the  proton  quasiparticles.  [92] 

Associated  with  the  centroid  variable  is  the  path  integral  centroid  density,  denoted  by 
Pc  (RH!,-,RHi^;i'^''^;R0i,-,R'O;vo)  present  case.  This  equilibrium  density,  which 

bears  the  closest  resemblance  to  the  classical  Boltzmann  density,  is  formally  calculated 
by  fixing  the  centroid  variables  in  Eq.  (2.12)  and  then  integrating  over  all  configurations 
of  the  discretized  quasiparticles,  subject  to  a  weighting  by  their  effective  Boltzmann-like 
configurational  factor.  While  the  centroid  picture  in  equilibrium  statistical  mechanics  is 
interesting  in  its  own  right  [93-97],  it  was  the  discovery  [88,89]  and  subsequent  proof  [90]  of 
the  classical-like  dynamical  properties  of  the  centroid  variable  that  have  made  the  present 
study  possible. 

The  basic  notion  of  CMD  is  that  the  position  or  velocity  correlation  function  for  quan¬ 
tum  particles  in  a  many-body  system  can  be  approximately  related  to  a  centroid  position 
or  velocity  correlation  function  obtained  by  running  classical-like  trajectories  on  the  effec¬ 
tive  quantum  centroid  potential  [88,96,89-91,97].  The  centroid  trajectories,  as  well  as  the 
trajectories  of  any  classical  degrees  of  freedom  such  as  the  oxygens,  are  generated  by  the 
classical-like  MD  equations 


(2.13) 

and 

mo  •  Ro(^)  =  (Rh\Ro)  > 

(2.14) 

where  the  terms  mo  and  mn  are  vectors  with  components  mo  and  mn. 

centroid  potential  is  formally  defined  as 

The  effective 

F(‘=)(rS\Ro)  =  -fcsTln  [pc  (R?,Ro)^  • 

(2.15) 

The  centroid  density  in  the  discrete  representation  is  simply  given  by 

p,  (<’.  Ro)  =/••■/  n  ^Rh  ^  (Rh’  -  Rt)  . 

i=l 

(2.16) 
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The  centroid  forces  (Fh^)  on  the  protonic  degrees  of  freedom  in  the  discretized  representation 


are  given  by 

(,)  _  /  •  •  •  J  ULi  dR^^6  -  Rgi)  V  H  exp  {-pVp} 

UL  -  R&)  exp{-/3Vp} 


(2.17) 


where,  from  the  cyclic  invariance  of  the  trace, 

V'H  =  ;^f;VR£,V'(Rl{>,Ro). 

^  i=l 

Similarly,  the  forces  for  the  classical  degrees  of  freedom  (Fq)  are  given  by 

J  ■  ■  ■  J  nr=i  rfRS*^  (Rg  -  Rh)  V'o  exp  {-/3Vp} 

/■■■lnf.idRil^S(R§’-R?,}exp{-0Vp} 

where 


(2.18) 


(2.19) 


V'o  =  ii;VR.,V(Rli>,R<,).  (2.20) 

^  i=l 

The  CMD  method  includes  the  effects  of  quantum  zero-point  energy  and  tunneling  in 
molecular  dynamics  simulations  directly  and  efficiently.  While  the  method  is  approximate, 
the  supporting  results  are  quite  compelling.  For  example,  CMD  has  proven  to  be  very 
accurate  in  several  simulations  [88,96,89-91,98]  and  it  has  been  justified  through  a  direct 
mathematical  analysis  [90],  as  well  as  from  a  variational  perspective  [88,96].  Furthermore, 
when  barriers  are  encountered  by  the  centroid  trajectories,  the  insights  from  the  earlier  PI— 
QTST  [59-62]  confirm  that  the  trajectories  will  overcome  such  barriers  with  the  accurate 
quantum  probability. 


III.  MODEL  PARAMETERIZATION 
A.  HsO^  dimer  parameterization 

An  important  issue,  of  course,  is  the  way  in  which  to  model  the  interactions  of  a  hydro- 
nium  with  a  given  water  molecule  in  the  first  solvation  shell.  To  this  end,  the  ground  state 
electronic  surface  on  which  the  proton  transfers  between  two  water  molecules  was  described 
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using  the  EVB  method  [53,54].  In  this  method,  two  distinct  diabatic  electronic  states  are 
used  (cf.  Fig.  3)  in  which  the  proton  is  bound  to  one  water  molecule  in  the  first  state  (Fu) 
and  to  the  other  water  molecule  in  the  second  state  (V22)-  Empirical  potentials  were  used 
for  the  stable  hydronium  states,  while  information  from  explicit  electronic  structure  calcu¬ 
lations  was  used  to  determine  the  off-diagonal  matrix  elements  which  allow  the  transitions 
between  the  diabatic  valence  bond  states. 

To  be  more  specific,  the  two  valence  bond  states  interact  through  an  off-diagonal  matrix 
element  Vi2  which  can  be  parameterized  as  a  function  of  a  subset  of  the  nuclear  coordinates. 
The  adiabatic  ground  state  potential  energy  function  (Kd)  foi'  the  system  nuclei  of  this 
simple  dimer  is  found  by  diagonalizing  the  2x2  EVB  matrix  at  each  configuration,  i.e., 

Vad  —  ^  (Vii  +  V22)  -  -\/ (Vii  -  ^22)^  +  (^•^) 

In  this  expression,  the  empirical  potential  for  a  given  diabatic  state  i  is  written  as 

v,i  =  v«t°*  +  v;K  +  Vp.,,  (3.2) 

where  Vpair  is  the  intermolecular  interaction  between  the  water  molecule  and  hydronium 
ion  and  (V^intra^)  is  the  intramolecular  potential  of  the  water  (hydronium).  The  in¬ 

tramolecular  potential  for  the  hydronium  was  composed  of  the  three  harmonic  OH  bonds 
(roHi)  and  three  HOH  angle-bend  (6>i)  terms  given  by, 

=  t  koH  {roH,  -rof+  h  (9,  -  Dof  .  (3.3) 

i=l 

The  equilibrium  bond  length  (ro),  bond  angle  (0o))  and  force  constants  {koH  and  kg)  were 
chosen  to  reproduce  the  four  distinct  vibrational  frequencies  of  the  hydronium  ion.  The  two 
low  frequency  bend  modes  and  V2  of  the  model  were  adjusted  to  reproduce  the  aqueous 
phase  values  of  1730  cm"^  and  1200  cm“^  [27].  The  two  high  frequency  symmetric  {ui)  and 
asymmetric  (1/3)  stretches  were  adjusted  to  match  the  experimental  gas  phase  values  of  3650 
cm'^  and  3730  cm"^  [99]  since  these  are  more  decoupled  from  the  condensed  phase  effects. 
The  resulting  normal  modes,  as  well  as  the  experimental  gas  phase  vibrational  frequencies 
are  listed  in  Table  I.  The  parameters  for  the  dimer  model  are  given  in  Table  II. 
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The  intramolecular  potential  for  the  “acceptor”  water  molecule  was  the  harmonic  pa¬ 
rameterization  of  Kutchiso  and  Morino  [100]  which  is  given  by  [101] 

^intra  =  E  [bHH  - 

i—1 

+  c  (boHi  +  boH2  -  SfeoFe,)  {b  HH  -  bHHeg) 

+  d  {bo Hi  -  boH,^)  {boH2  -  boHe^  .  (3-4) 

where  612  in  the  bond  length  between  atoms  1  and  2.  The  equilibrium  bond  lengths  are  given 
in  Table  11.  The  intermolecular  interactions  between  the  water  molecule  and  the  hydronium 
ion  were  taken  to  be  pairwise  Lennard-Jones  (LJ)  and  Coulomb  interactions  between  the 
various  ionic  sites.  These  interactions  may  be  written  as 


where  e  and  cr  are  the  usual  LJ  parameters  and  q-m  {Qn)  are  the  partial  charges  on  the  mth 
(nth)  atom  of  the  hydronium  (water).  The  LJ  interactions  are  only  between  the  oxygen 
atoms.  The  partial  charges  were  taken  from  the  Mulliken  charge  population  analysis  in  the 
gas  phase  density  functional  electronic  structure  calculation  of  Wei  and  Salahub  [102]  for  the 
H5OJ  complex.  The  partial  charges,  as  well  as  the  LJ  parameters,  are  listed  in  Table  III. 

As  stated  previously,  the  barrier  in  HsO^  is  low  to  nonexistent  in  the  gas  phase,  de¬ 
pending  on  the  0-0  distance.  The  off-diagonal  term  which  couples  the  diabatic  states  was 
parameterized  to,  in  part,  reproduce  this  barrier.  The  parametric  form  for  V12  was  chosen 
to  be  given  by 

V12  =  V"i2  exp  -aoo  {Roo  -  -Roo)  > 

where  is  65.89  kcal  mol~\  aoo  is  8  and  is  2-6  A.  This  parameterization  yields 
a  Vad  in  Eq.  (3.1)  that  reproduces  reasonably  well  the  gas  phase  ab  initio  potential  energy 
surface  of  Tortonda,  et  al.  [36].  The  potential  surface  of  the  latter  study  was  generated  by 
calculating  the  ab  initio  potential  energy  at  constrained  00  and  OH  distances  while  allowing 
for  relaxation  of  all  the  other  coordinates  of  the  complex  to  their  potential  energy  minimum. 


The  analogous  surface  was  generated  for  the  empirical  potential  energy  function  in  Eq.  (3.1) 
through  a  MC  annealing  calculation  while  constraining  the  same  coordinates.  The  resulting 
potential  surface  is  shown  in  Figure  4.  The  relative  separation  of  the  minima  and  height  of 
the  transition  state  are  in  good  agreement  with  results  of  the  electronic  structure  calculation. 
The  height  of  the  barrier  is  approximately  0.35  kcal  mol"’'  which  compares  well  with  the 
barrier  of  0.3  kcal  mole"^  calculated  by  Tortonda,  et  al.  [36],  though  the  equilibrium  0-0 
separation  is  somewhat  larger  in  agreement  with  the  ab  initio  MD  result  of  Tuckerman,  et 
al.  [45] 


B.  Solvent-  H5O2  Complex  Interactions 

The  usual  approach  in  the  EVB  method  is  to  add  the  complex-solvent  interactions 
directly  to  each  diabatic  state  Vii  and  V22  iii  Eq.  (3.2).  The  potential  for  the  entire  solvent 
+  solute  system  is  then  found  by  diagonalizing  the  2x2  diabatic  matrix.  An  alternative 
approach  is  taken  here.  The  gas  phase  adiabatic  surface  of  Eq.  3.1  is  assumed  to  be  a 
good  description  of  the  intramolecular  degrees  of  freedom  of  the  complex  and  relatively 
unaffected  by  the  solvent.  Hence  the  adiabatic  state  of  the  solute  is  a  function  of  the  solute 
coordinates  alone  and  does  not  depend  on  the  solvent.  The  total  potential  energy  including 
the  solvent-solvent  and  solvent-complex  interactions  is  thereby  written  as 

Id  =  Fad  +  (R-h,  R-o)  +  (Rh,  Rq,  r) ,  (3.7) 

where  V55  contains  the  solvent-solvent  interactions  as  well  as  the  intramolecular  potential 
energy  of  the  solvent  molecules.  The  term  Vsc  is  the  potential  energy  of  the  solvent-complex 
coupling.  For  the  water  solvent,  the  SPC  water  model  [103]  with  the  intramolecular  harmonic 
parameterization  of  Kutchiso  and  Morino  [100,101]  was  used  [cf.  Eq.  (3.4)].  The  equilibrium 
and  dynamical  properties  of  the  quantized  SPC/F  water  model  will  be  fully  explored  in  a 
future  publication  [104]. 

As  the  excess  proton  transfers  between  two  molecular  water  hosts,  there  is  considerable 
redistribution  of  charge  as  the  net  unit  positive  charge  shifts  its  distribution  over  the  donor 
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hydronium  to  the  acceptor  water.  In  order  to  mimic  this  charge  redistribution,  the  electro¬ 
static  atomic  site  charges  of  the  complex  were  parameterized  as  functions  of  the  reaction 
coordinate,  the  asymmetric  stretch  of  the  transferring  proton  [cf.  Eq.  (2.4)].  To  this  end,  the 
ab  initio  study  of  Wei  and  Salahub  [102]  was  used  for  parameterizing  the  partial  charges  of 
the  H5O2  complex  with  which  the  solvent  interacts.  In  the  latter  study,  a  Mulliken  charge 
population  analysis  was  carried  out  for  all  the  atoms  of  the  H5O2  complex  and  charges 
were  calculated  at  different  positions  along  the  proton  transfer  reaction  coordinate.  The 
dependence  of  these  charges  on  the  reaction  coordinate  was  thus  reproduced  in  the  empiri¬ 
cal  model.  The  following  equations  were  used  to  reproduce  the  charge  variation.  Using  the 
atomic  site  labeling  in  Fig.  2,  the  partial  charges  for  the  protons  of  the  complex  were  given 
by 


^Hi  =  /  {Qasym)  +  [1  /  {Qasym)]  ? 

(3.8) 

^H2  “  /  {Qasym)  +  [l  “  /  {Qasym)]  ? 

(3.9) 

—  9  {Qasym)  ? 

(3.10) 

eH4  =  {Qasym)  +  [l  “*  /  {Qasym)]  i 

(3.11) 

{qasym)  +  [l  ”  /  (^asym)]  • 

(3.12) 

The  Coulomb  charges  of  the  oxygens  were  given  by 

eoi  =  {Qasym)  +  [l  "  /  {lasym)]  +  ^  “  9  {Qasym)]  ,  (3.13) 

e02  =  eo'°/  {Qasym)  +  60'*^^  /  iQasym)]  +  ^  [e®  “  9  {Qasym)]  ■  (3-14) 

The  two  charge  switching  functions  /  {qasym)  and  g  {qasym)  were  given  by 

/  {Qasym)  =  \  ~  {Qasym/ Qsw)]  ,  (3-15) 

and 

9  {Qasym)  ^swQasym  Qmin'  (3.16) 
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The  parameters  used  in  the  above  equations  are  given  in  Table  IV.  The  above  empiri¬ 
cal  equations  and  ah  initio  charges  of  Wei  and  Salahub  are  plotted  as  a  function  of  the 
asymmetric  stretch  coordinate  in  Fig.  5. 


C.  Simulation  Details 

The  discretized  partition  function  Qp  in  Eq.  (2.9)  converges  rigorously  to  the  quantum 
limit  as  P  — >  00.  However,  in  practice  a  finite  value  of  P  is  used  so  that  the  thermodynamic 
properties  of  interest  are  sufl&ciently  converged  to  their  quantum  value.  In  the  present  case,  a 
MC  study  was  performed  to  select  a  suitable  value  of  P  [104,105].  This  study  followed  along 
the  lines  of  one  suggested  by  Kuharski  and  Rossky  [106]:  It  consisted  of  two  solvent  water 
molecules  which  interacted  using  the  harmonic  SPC/F  potential  (described  in  Sec.  IIIB) 
and  an  additional  quadratic  potential  between  their  centers  of  mass.  This  potential  had  a 
minimiiTn  at  2.85  A  and  a  frequency  u;  =  26  ps“\  creating  an  environment  similar  to  the 
bulk  solvent  and  typical  of  water  librational  motion.  A  total  of  10®  trial  MC  moves  were 
carried  out  and  a  number  of  different  properties  were  examined  as  a  function  of  P  including 
radial  distribution  functions  and  the  average  intramolecular  energy.  The  value  P  —  25  was 
found  to  yield  sufficient  convergence  for  these  properties. 

The  fictitious  mass  for  the  PIMD  algorithm  based  on  Eq.  (2.11)  should  be  chosen 
to  be  as  small  as  possible  in  order  to  effect  rapid  sampling.  However,  if  it  is  too  small, 
ergodicity  problems  will  occur  and/or  numerical  integration  of  the  equations  of  motion  will 
be  difficult  due  to  the  small  time  step  necessary  to  integrate  the  high  frequency  motion  of  the 
proton  quasiparticles.  As  a  result,  the  fictitious  mass  was  chosen  by  a  normal  mode  analysis 
of  a  single  quantized  water  molecule  with  the  harmonic  intramolecular  SPC/F  potential  and 
P  =  25.  The  mass  was  set  to  200  au  so  that  the  highest  frequency  mode  was  below  6000 
cm“h 

The  Velocity  Verlet  MD  [107]  method  with  an  implicit  iterative  algorithm  described  by 
Tuckerman,  et  al.  [108,109]  was  used  to  integrate  the  PIMD  Hamiltonian.  Two  Nose-Hoover 
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oscillators  [110-112]  were  used  for  equilibration,  one  “attached”  to  the  oxygens,  the  other 
to  the  hydrogen  quasiparticles.  It  was  verified  that  these  two  Nose  oscillators  were  sufiicient 
to  overconae  ergodicity  problems  associated  with  PIMD  calculations  [113].  The  Nose  masses 
{Qr]o  the  oxygens  and  Qrj^  for  the  protons)  were  chosen  to  be  [114] 


Qvo 


Qrju  — 


SNo-3 

^wat^ 

SPiVn 

:/? 


hydr 


(3.17a) 

(3.17b) 


where  and  cohyd  are  the  characteristic  frequencies  of  the  system  and  were  set  to  a;TOat=680 
cm”^  and  Whi/d^SGOO  cm“^  if  the  protons  were  treated  classically  or  uJhyd  = 
if  the  protons  were  treated  quantum  mechanically.  The  time  step  used  for  both  the  classical 
and  PIMD  simulations  was  0.5  fs. 


D.  Centroid  Molecular  Dynamics  Algorithm 

A  number  of  algorithms  for  solving  the  CMD  equations  [cf.  Eqs.  (2.13)-(2.17)]  have  been 
developed  [91].  This  task  is  not  entirely  trivial  since  the  centroid  potential  y(‘^^(RH\Ro) 
is  actually  a  quantum  potential  of  mean  force,  requiring  path  integral  averaging  to  find  the 
centroid  force  at  each  time  step  [cf.  Eqs.  (2.16)-(2.17)].  In  the  approach  taken  here  the 
“natural”  CMD  time  step  was  broken  into  Nmd  smaller  time  steps.  At  each  of  these  small 
time  steps  a  PIMD  calculation  was  run  to  obtain  the  centroid  force.  The  centroid  forces 
are  the  forces  felt  by  the  centroids  of  the  quantum  degrees  of  freedom  and  the  positions 
of  the  classical  degrees  of  freedom.  Hence,  the  PIMD  calculation  was  run  with  constraints 
on  the  Cartesian  positions  of  the  centroids  [115]  without  moving  the  classical  degrees  of 
freedom.  The  forces  were  averaged  over  a  reasonable  number  of  PIMD  time  steps  followed 
by  a  MD  move  of  the  proton  centroids  and  the  classical  oxygen  coordinates.  This  “on  the 
fly”  algorithm  provides  a  feasible  alternative  to  a  brute  force  approach  for  obtaining  the 
centroid  force.  The  CMD  time  step  used  was  0.05  fs  and  the  PIMD  time  step  used  was  0.5 
fs.  The  centroid  force  was  averaged  over  five  PIMD  time  steps  (e.g.,  five  PIMD  timesteps 
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were  used  for  every  CMD  time  step).  The  simulations  were  run  at  300  K  with  125  molecules 
(123  water  molecules  and  the  HgOj  complex)  at  the  density  of  water  (1.00  g/mL).  The 
interactions  were  tapered  at  half  the  box  length  with  a  smooth  spherical  cutoff. 

IV.  RESULTS 

A.  Radial  Distribution  Functions 

The  microscopic  structure  of  the  solvent  around  the  HsO^  complex  was  first  analyzed 
using  radial  distribution  functions.  In  Fig.  6,  the  radial  distribution  of  oxygens  about  each 
oxygen  of  the  complex  shows  a  sharp  maximum  at  2.5  A.  The  positioning  of  the  peak  is  in 
agreement  with  the  experimental  thermal  neutron  and  X-ray  scattering  studies  of  Triolo,  et 
al.  [28].  If  one  regards  the  H3O+  as  a  chemically  distinct  species,  the  maximum  corresponds 
to  the  three  water  molecules  that  hydrogen  bond  to  the  three  hydrogens  of  the  HsO'*'  in  the 
HgOi"  complex.  The  peak  is  moved  in  and  sharper  than  the  oxygen-oxygen  peak  of  liquid 
water  which  is  shown  for  comparison.  However,  the  peak  is  split  in  this  case  because  one 
of  the  three  oxygens  is  not  a  water  oxygen  but  a  companion  oxygen  in  the  H5O2  complex 
which  is  somewhat  closer  than  the  oxygens  of  the  two  water  molecules  that  hydrogen  bond 
to  a  given  oxygen  in  the  HsO^  complex.  This  splitting  was  also  seen  in  a,b  initio  simulations 
of  Tuckerman,  et  al.  [31].  In  the  latter  simulations,  two  different  radial  distributions  were 
calculated.  In  the  first  simulation  no  splitting  was  observed  which  corresponded  to  the 
symmetrical  HgOj  structure  in  water.  In  the  second  distribution  function  the  peak  was 
split,  indicating  that  one  of  the  water  molecules  was  moved  in  and  the  H5O2  complex  was 
present.  Since  the  distribution  functions  are  averaged  over  very  long  time  intervals  in  the 
present  study  the  first  peak  is  split  as  the  two  structures  interchange.  It  should  be  noted 
that  the  charge  switching  parameterization  of  the  empirical  model  allows  the  acceptor  water 
to  have  a  charge  distribution  which  is  almost  the  same  as  that  of  SPC/F  water  when  the 
transferring  proton  is  close  to  the  donor  hydronium.  In  addition,  the  model  is  parameterized 
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in  such  a  way  that  the  charge  distribution  of  the  donor  hydronium  is  symmetrical  over  the 
three  atoms  of  the  H3O+  .  These  features  allow  for  the  dynamical  interchange  of  the 
EgOt  structure  and  the  HgO^  complex.  Both  of  these  were  in  fact  observed  in  the  present 
simulations  using  computer  visualization  packages.  The  splitting  is  enhanced  in  the  quantum 
case  because  the  average  oxygen-oxygen  distance  is  greater  for  the  two  water  molecules  that 
hydrogen  bond  to  the  oxygen  of  the  complex.  This  is  due  to  the  fact  the  average  OH  bond 
length  of  the  four  hydrogens  that  do  not  participate  in  proton  transfer  in  the  HsO^  complex 
is  greater  in  the  quantum  case.  The  OH  bond  length  is  1.048  ±  0.002  A  in  the  classical  and 
1.068  ±  0.001  A  in  the  quantum  case.  The  average  00  distance  of  the  oxygens  in  the 
H5OJ  complex  is  2.4862  ±  0.0007  in  the  classical  and  2.4860  ±  0.0006  in  the  quantum  cases. 
Both  numbers  are  in  good  agreement  with  the  experimentally  observed  distance  of  2.5  A  [28]. 
The  coordination  number  at  the  first  minimum  of  the  oxygen-oxygen  radial  distribution  is 
3.3,  while  the  experimentally  determined  coordination  number  is  four  [28].  It  should  be 
noted,  however,  that  the  coordination  number  for  the  oxygen  in  the  H3O+  at  the  minimum 
of  the  oxygen— oxygen  radial  distribution  for  liquid  water  (~  3.3  A)  is  four.  Therefore,  the 
coordination  number  of  the  H30'^  when  evaluated  at  the  first  solvation  shell  radius  of  liquid 
water  is  more  in  accord  with  experiment. 

Further  insight  into  the  microscopic  structure  of  the  solvent  can  be  gained  by  examining 
the  distribution  of  hydrogens  with  respect  to  the  oxygens  of  the  complex  (cf.  Fig.  7).  The 
transferring  proton  has  been  included  in  this  distribution  function,  but  the  four  hydrogens 
covalently  bound  to  the  complex  have  not.  In  the  classical  treatment  of  the  protons,  the 
peak  due  to  the  transferring  proton  is  distributed  between  1.0  and  1.5  A  and  is  split  due 
to  the  fact  that  there  is  a  small  effective  barrier  (double  well  potential)  along  the  proton 
transfer  coordinate.  In  the  quantum  case,  the  peak  is  not  split  due  to  the  fact  that  the 
proton  tunnels  and  has  a  large  zero  point  energy  above  the  classical  gas  phase  barrier  (cf. 
the  following  section).  The  peak  at  1.9  A  present  in  the  water  OH  distribution  function 
is  not  present  in  the  hydronium  case,  indicating  that  the  protons  of  the  waters  in  the  first 
solvation  shell  are  oriented  outwards  with  respect  to  the  complex.  The  classical  distribution 
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function  calculated  here  is  quantitatively  the  same  as  that  found  by  Tuckerman,  et  ah  [31] 
using  ab  initio  molecular  dynamics,  although  the  statistical  error  is  much  greater  in  the 
ah  initio  simulation.  This  confirms  that  the  parameterization  of  the  empirical  model  is 
quite  successful  when  compared  with  more  sophisticated,  but  much  more  costly,  dynamical 
electronic  structure  methods.  The  present  result  also  indicates  that  the  quantum  effects  for 
the  transferring  proton  are  significant,  qualitatively  changing  the  distribution  function. 

The  distribution  function  of  the  transferring  proton  with  the  protons  of  the  complex  and 
protons  of  the  nearby  water  molecules  shows  a  large  maximum  corresponding  to  the  four 
protons  of  the  complex  (cf.  Fig.  8).  Again,  the  peak  is  split  in  the  classical  case  due  to 
the  barrier  along  the  proton  transfer  coordinate.  The  distribution  shows  enhanced  ordering 
of  the  solvent  at  longer  distances  due  to  the  overall  charge  of  the  complex  which  aligns 
the  solvent  water  molecules.  All  of  the  radial  distribution  functions  discussed  above  were 
averaged  over  600  ps. 


B.  Asymmetric  Stretch  Trajectories 

Representative  trajectories  of  the  proton  transfer  asymmetric  stretch  coordinate  are 
shown  in  Fig.  9,  wher  the  classical  (Fig.  9a)  and  centroid  (Fig.  9b)  trajectories  of  the 
asymmetric  stretch  coordinate  are  plotted  as  functions  of  time.  The  high  frequency  oscilla¬ 
tions  of  the  quantum  centroid  trajectory  are  of  larger  amplitude  than  those  of  the  classical 
variable.  In  addition,  the  classical  transfer  coordinate  sometimes  resides  on  either  side  of 
the  classical  transition  state  {qasym  —  0)  whereas  the  centroid  trajectory  samples  the  entire 
range  of  qasym  in  smooth  fashion,  indicating  a  lower  or  non-existent  activation  free  energy 
along  the  proton  transfer  coordinate  in  the  quantum  case.  In  the  classical  case,  the  value  of 
qasym  corresponds  to  the  case  of  the  H9O4  structure  about  sixty  percent  of  the  time  and  to 
the  H5O2  complex,  in  which  the  proton  resonates  between  water  molecules,  around  forty 
percent  of  the  time.  This  classical  result  is  again  in  qualitative  agreement  with  the  ab  initio 
MD  results  of  Tuckerman,  et  al.  [31]  for  classical  nuclear  dynamics,  but  the  effects  of  nuclear 
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quantization  are  again  seen  to  be  large  in  the  present  study. 


C.  Quantum  Dynamics  of  Proton  Migration 

The  activation  energy  of  proton  transfer  in  the  HsO^  complex  was  calculated  in  both 
the  classical  TST  (i.e.  all  the  protons  in  the  system  were  classical)  and  the  quantum  PI- 
QTST.  As  can  be  seen  from  Fig.  10,  in  the  classical  limit  the  activation  energy  is  about 
0.51  kcal  mol~^,  which  is  approximately  0.9  fc^T,  and  arises  from  the  solvent  orientational 
polarization  reorganization.  There  is  no  centroid  activation  free  energy  in  the  quantum  case. 
The  dependence  of  the  activation  free  energy  on  system  size  was  also  checked;  curves  for 
both  123  and  510  water  molecules  are  shown  in  Fig.  10.  Proton  “sharing”  between  two  water 
molecules  in  the  quantized  H5OJ  complex  is  clearly  an  activationless  or  nearly  activationless 
process  and  certainly  not  the  rate  limiting  step  to  proton  migration  in  water.  This  result 
could  have  been  predicted  from  the  character  of  the  centroid  trajectory  in  Fig.  10b. 

Interestingly,  the  activation  energy  of  the  rate  limiting  step  for  proton  transfer  in  water 
has  been  measured  in  experiments  by  Luz  and  Meiboom  [8]  to  be  approximately  2.4 
kcal  mol~^.  Furthermore,  the  mean  residence  time  of  a  proton  with  a  given  water  molecule 
is  around  1.5  ps  [8,9,13].  The  inverse  of  this  number  is  the  rate  of  proton  transfer  and  is 
0.69  ps“^  [8].  The  rate  is  clearly  much  smaller  than  the  rate  of  proton  transfer  through 
the  strong  hydrogen  bond  in  the  H5OJ  complex  as  evidenced  by  Fig.  9b.  Extended  proton 
migration  must  therefore  occur  in  a  concerted  fashion,  perhaps  requiring  both  the  correct 
solvent  fluctuations  in  the  second  solvation  shell  of  the  H5OJ  complex  and  an  inward 
fluctuation  of  the  oxygen-oxygen  distance  between  the  complex  and  another  water  molecule 
(i.e.,  the  coordinate  Rc  in  Fig.  2).  The  inward  fluctuation  of  this  oxygen— oxygen  distance 
creates  a  “special”  bond  wherein  the  resonating  proton  “flips”  its  identity  (i.e.,  from  proton 
three  to  proton  four  in  Fig.  2),  thus  forming  a  new  H5O2  complex  which  includes  the 
inward  fluctuating  water  molecule.  It  is  this  switching  of  the  “special”  resonating  proton 
hydrogen  bond  between  pairs  of  oxygen  atoms  that  could  give  rise  to  the  Gr6tthus-type 
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proton  migration  mechanism.  [45] 

In  order  to  gain  better  insight  into  this  process,  the  rate  of  inward  00  fluctuations  of  the 
four  waters  closest  to  the  HsO^  complex  was  calculated.  Xhis  rate  was  also  correlated  with 
the  probability  that  the  four  water  molecules  were  sufficiently  hydrogen  bonded  to  other 
water  molecules.  Both  conditions  are  necessary  for  a  flip  of  the  “special”  bond  to  occur 
and  thus  for  extended  proton  transfer  to  occur.  In  order  to  define  the  energy  criterion  of  a 
hydrogen  bond  between  a  water  in  the  first  solvation  shell  of  the  H5OJ  complex  and  the 
nearby  solvent  molecules,  the  distribution  of  pairwise  energies  was  calculated  for  the  four 
waters  closest  to  the  H5O2  complex  with  the  other  water  molecules  of  the  solvent.  The 
classical  and  quantum  distributions  are  shown  in  Fig.  11.  The  hump  on  the  shoulder  of  the 
distribution  is  due  to  hydrogen  bonding. 

The  rate  of  inward  oxygen-oxygen  fluctuations  was  calculated  by  first  “tagging”  the 
four  water  molecules  closest  to  the  oxygen-oxygen  complex.  A  dividing  surface  was  placed 
at  a  distance  Rc  and  the  rate  at  which  the  oxygen-oxygen  distance  crossed  the  dividing 
surface  was  measured.  A  crossing  of  the  dividing  of  the  surface  was  only  counted  if  the 
water  molecule  had  two  pairwise  interactions  with  the  solvent  molecules  that  were  below 
—4.0  kcal  mol“h  This  definition  of  a  hydrogen  bond  is  reasonable  given  the  distribution 
in  Fig.  11  and  that  the  average  energy  of  a  hydrogen  bond  for  classical  rigid  classical  SPC 
water  [116]  is  -4.34  kcal  mole“h  The  resulting  rate  should  represent  an  upper  bound  on  the 
rate  of  flipping  of  the  “special”  bond  and  hence  the  rate  of  proton  transfer.  The  resulting 
rates  are  plotted  as  a  function  of  the  oxygen-oxygen  distance  Rc  in  Figure  12.  As  can  be 
seen  from  that  figure,  the  overall  rate  is  significantly  lower  in  the  quantum  case  (which  was 
calculated  using  CMD).  The  quantum  rate  in  the  range  of  Rc  =  2.5  —  2.6  A  is  in  good 
agreement  with  the  experimentally  measured  room  temperature  rate  of  the  proton  hop  of 
0.69  ps~^  [8].  This  distance  range  is  considered  to  be  typical  of  the  oxygen-oxygen  distances 
at  which  the  “flip”  of  the  special  complex  H-bond  could  take  place.  The  classical  rate  in 
the  same  range  of  Rc  is  approximately  10-20  ps~h  The  quantum  rate  calculation  shown 
in  Figure  12  was  averaged  over  20  ps  of  a  CMD  simulation,  while  the  classical  rate  was 
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averaged  over  40  ps.  Clearly,  the  quantum  effects  on  the  nuclear  motion  are  significant  in 
this  case. 


D.  Zundel  Polarization  and  the  Continuum 

Strong  hydrogen  bonds  are  characterized  by  low  to  nonexistent  barriers  for  proton  trans¬ 
fer,  enthalpies  of  9-15  kcal  mor\  and  short  equilibrium  distances  (e.g.  0-0  distances  less 
then  2.6  A).  The  spectral  signature  of  strong  H-bonds  in  solution  (and  acid  crystals)  is  a 
strong  absorption  peaked  at  1500-2000  cm“^  and  extending  over  a  range  of  1000-1500  cm"k 
This  absorption  is  known  as  the  continuum  and  is  observed  in  many  liquid  systems  (not  just 
water)  in  which  symmetric  H— bonds  of  the  type  BH’*’  •  •  -B  or  (AH-  •  -A)  are  present.  [29] 
This  phenomenon  is  a  protonic  rather  than  an  electronic  polarization  effect.  The  polarization 
of  a  proton  in  a  strong  hydrogen  has  been  estimated  to  be  one  to  two  orders  of  magnitude 
larger  than  electronic  polarizabilities  of  molecules  [38-40].  Two  mechanisms  are  thought  to 
be  responsible  for  the  continuum: 

1.  Zundel  polarization  which  is  also  called  the  direct  mechanism.  Strong  electrostatic 
coupling  of  the  solvent  gives  rise  to  deformations  of  the  proton  transfer  potential. 
There  is  a  large  distribution  of  fields  which  gives  rise  to  a  continuity  of  energy  level 
differences.  In  addition,  there  is  a  large  distribution  of  equilibrium  H  bond  lengths, 
causing  additional  “smearing”  of  the  continuum  [29]. 

2.  The  indirect  mechanism,  which  is  the  coupling  of  the  OH  vibration  to  the  low  frequency 
anharmonic  0-0  vibration  of  the  proton  transfer  complex  (or  whichever  two  atoms 
are  sharing  the  proton  in  a  strong  hydrogen  bond  [117-121]). 

A  discussion  of  the  effects  of  a  static  external  field  on  the  absorption  spectrum  of  the 
H5OJ  complex  is  given  by  Janoschek  [122].  One  of  the  conclusions  of  this  work  and  re¬ 
lated  treatises  [38-40]  is  that  the  electric  field  produced  by  the  surrounding  polar  solvent  is 
sufficient  to  lead  to  considerable  broadening  of  the  IR  absorption  spectrum  in  strongly  hy- 
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drogen  bonded  (i.e.  low  barrier)  systems.  A  rigorous  study  to  obtain  the  infrared  absorption 
spectrum  would  require  finding  the  line  shape  function  I{uj)  which  is  given  by  [123] 

/(cj)  =  r  dte^^  {MKt))  ,  (4-1) 

J — oo 

where  (/r(0)/r(t))  is  the  macroscopic  quantum  dipole  autocorrelation  function.  Since  only  a 
qualitative  comparison  is  made  here,  the  quantupi  velocity  autocorrelation  function  of  the 
O1H3  bond  is  used  instead,  i.e., 

I,{u)  =  r  {vo,Hs{0)voMt))  ■  (4-2) 

J-oo 

Since  this  vibrational  mode  is  directly  involved  in  the  transfer  of  the  proton,  its  power 
spectrum  should  contain  most  of  the  qualitative  features  of  the  continuum  IR  spectrum  of  a 
strong  acid  dissolved  in  water.  Within  the  framework  of  the  theory  of  CMD,  the  quantum 
correlation  function  is  obtained  from  the  centroid  correlation  function  (which  is  calculated 
from  computer  simulation)  using  the  following  expression  [89,90] 

Iv{u))  =  {hl3uj/2)  [coth  {h(3uj/2)  +  1]  (4-3) 

where  is  the  Fourier  transform  of  the  centroid  velocity  correlation  function.  The 

power  spectrum  is  plotted  in  Figure  13  and  was  obtained  by  averaging  over  four  ps  of  a 
CMD  run.  The  strong  absorption  between  1000  and  2500  cm“^  is  typical  of  strong  acids  in 
solution.  Some  residual  absorption  at  3600  cm“^  is  present  from  the  OH  covalent  stretch.  As 
stated  earlier  the  absorption  is  thought  to  occur  via  two  mechanisms,  the  direct  mechanism 
due  to  solvent  deformations  of  the  proton  potential  giving  rise  to  a  continuum  of  transition 
and  the  indirect  mechanism  due  to  coupling  of  the  OH  stretch  to  the  low  frequency  0-0 
vibration.  Differentiation  between  the  two  mechanisms  could  be  made  by  constraining  the 
oxygen-oxygen  distance  in  the  complex,  thus  eliminating  the  contribution  of  the  indirect 
mechanism.  Future  work  will  focus  on  this.  The  spectrum  presented  here  is  also  in  good 
agreement  with  simulations  of  condensed  phase  proton  transfer  in  a  model  system  by  Borgis 
et  al.  [124]. 
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V.  CONCLUDING  REMARKS 


In  this  paper,  structural  and  dynamical  studies  of  an  excess  proton  in  water  have  been 
carried  out.  The  coordination  number  of  the  H3O+  species  was  found  to  be  in  good  agree¬ 
ment  with  experiment  and  the  classical  radial  distribution  functions  of  the  HsO^  species  are 
in  accord  with  those  found  from  ab  initio  MD  methods.  Importantly,  the  nuclear  quantum 
effects  are  found  to  be  large,  accordingly  the  proton  transfer  along  the  strong  hydrogen  bond 
in  the  H5OJ  complex  is  found  to  be  a  quantum  activationless  process.  The  strong  hydrogen 
bond  in  the  complex  is  considered  to  be  a  “special”  bond  [31]  characterized  by  a  smaller 
than  normal  oxygen-oxygen  distance  (~  2.5  A)  which  can  lead  to  the  rapid  interchange  of 
the  chemical  and  hydrogen  bonds  as  the  proton  moves  back  and  forth.  As  a  result  of  our 
dynamical  study  of  the  conditions  necessary  for  this  to  occur  ,  the  activation  energy  and  rate 
constant  measured  by  Luz  and  Meiboom  [8]  would  seem  to  correspond  to  the  flipping  of  the 
special  bond  to  another  pair  of  oxygens.  The  flipping  of  the  special  bond  corresponds  to  the 
simplest  example  of  the  Grotthus  mechanism  of  proton  transfer,  and  it  depends  on  the  dy¬ 
namics  of  the  water  molecules  in  the  second  solvation  shell,  of  the  proton  transfer  complex, 
as  well  as  the  oxygen-oxygen  distance  fluctuations  of  water  molecules  that  hydrogen  bond 
to  the  complex.  Analysis  of  the  quantum  dynamics  of  the  water  molecules  which  hydrogen 
bond  to  the  HsO^  complex  show  that  the  rate  of  inward  oxygen-oxygen  fluctuations  of 
these  molecules,  when  correlated  with  solvent  hydrogen  bonding  from  water  molecules  in 
the  second  solvation  shell,  is  in  agreement  with  the  experimentally  measured  rate  of  proton 
transfer.  The  nuclear  quantum  effects  in  this  process  are  significant. 

Finally,  the  quantum  power  spectrum  of  the  proton  transfer  mode  of  the  HsO^  complex 
in  water  was  found  to  be  in  qualitative  agreement  with  the  the  IR  spectrum  of  strong 
hydrogen  bonding  systems,  showing  a  characteristic  broad  absorption.  Future  work  will 
be  devoted  to  exploration  of  whether  the  direct  or  indirect  mechanism  predominates  in 
the  calculated  spectrum,  as  well  as  to  actual  quantum  simulations  of  the  extended  proton 
migration  process. 
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TABLES 


TABLE  1.  Experimental  [27,99]  and  model  [125]  vibrational  frequencies  of  the  hydronium  ion 
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TABLE  II.  Intramolecular  water  and  hydronium  parameters. 


Parameter 

Value 

koH 

546  kcal  mol“^ 

ro 

0.98  A 

ke 

44.55  kcal  mol“^  radian"^ 

Oo 

110  degrees 

pvj 

2.566  A 

Dyj 

0.708  mdyn  A 

bOHeq 

1.0  A 

HOH  angle 

109.47  degrees 

b 

2.283  mdyn  A“^ 

c 

-1.469  mdyn  A“^. 

d 

0.776  mdyn  A“^ 
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TABLE  III.  Intermolecular  potential  parameters  for  the  hydronium  and  water. 


0-0 

H2O 

H3O+ 

e/kB  (K) 

c{k) 

90  (e) 

9H  (e) 

90  (e) 

9H  (e) 

78.22 

3.165 

-0.83826 

0.45212 

-0.641 

0.52518 
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TABLE  IV.  Charge  switching  parameters  used  in  Eg.  (3.12)-(3.16). 


Parameter 

Value 

pH30+ 

-0.64151  e 

pH20 

-0.83826  e 

H3O+ 

0.52804  e 

.H2O 

0.45212  e 

e-m 

0.51945  e 

^min 

0.58575  e 

Otsw 

10.8396  e 

Qsw 

0.3175  A 
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FIGURES 

FIG.  1.  Illustration  of  the  Grotthus  mechanism.  Charge  transfer  takes  place  via  an  interchange 
of  chemical  and  hydrogen  bonds. 

FIG.  2.  Schematic  of  the  proton  transfer  reaction  coordinate  and  atoms  of  the  proton  transfer 
complex.  The  transferring  proton  is  labeled  H3  and  ri  (r2)  is  the  distance  between  H3  and  Oi 
(O2).  Also  illustrated  are  four  water  molecules  adjacent  to  the  complex  and  the  critical  second 
0-0  (i?c)  distance  (cf.  Sec.  IV  C) . 

FIG.  3.  Diabatic  states  of  the  HsO^  complex. 

FIG.  4.  Gas  phase  empirical  potential  energy  surface  of  the  H5O2  complex  as  a  function  of 
the  Oi— O2  distance  and  Oi— H3  distance.  The  contours  are  in  1.0  kcal  mol  ^  increments. 

FIG.  5.  Plot  of  the  ab  initio  gas  phase  charges  of  atoms  of  the  H5O2  complex  and  empirical 
charge  switching  functions  used  to  reproduce  the  ab  initio  results.  Charges  are  plotted  as  a  function 
of  the  proton  transfer  reaction  coordinate.  The  charge  on  H2  is  the  same  as  Hi.  Charges  on  H4 
and  H5  are  the  mirror  image  of  Hi,  reflected  about  qasym  =  0.  Similarly,  the  charge  on  O2  is  the 
mirror  image  of  Oi. 

FIG.  6.  Radial  distributions  of  oxygen  atoms  of  the  H5O  J  complex  with  respect  to  the  solvent 
and  complex  oxygen  atoms.  Radial  distributions  are  shown  for  the  quantum  (solid)  and  classical 
(dashed  line)  cases.  The  classical  oxygen-oxygen  distribution  function  (dot-dashed)  of  SPC/F 
water  is  also  shown  for  comparison. 

FIG.  7.  Radial  distribution  of  water  hydrogens  as  well  as  the  transferring  proton  with  respect 
to  the  oxygen  atoms  of  the  complex  (quantum-solid,  classical-dashed).  The  classical  SPC/F  oxy¬ 
gen-hydrogen  distribution  function  of  liquid  water  is  also  shown  for  comparison. 
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FIG.  8.  Radial  distribution  of  water  hydrogens  and  covalently  bound  complex  hydrogens  with 
respect  to  the  transferring  proton  of  the  complex.  The  quantum  distribution  is  the  solid  line  while 
the  classical  the  dashed  line.  The  classical  hydrogen-hydrogen  distribution  function  of  liquid  water 
(SPC/F)  is  also  shown  for  comparison. 

FIG.  9.  Trajectory  of  the  transferring  proton  along  the  asymmetric  stretch  coordinate  m  the 
(9a)  classical  and  (9b)  quantum  centroid  cases. 

FIG.  10.  Activation  free  energy  curves  as  a  function  of  the  proton  asymmetric  stretch  coordi¬ 
nate.  The  classical  activation  energy  curve  is  shown  in  10a,  while  the  quantum  activation  energy 
curve  is  shown  in  10b.  Both  free  energy  curves  have  been  calculated  with  system  sizes  of  510  and 
123  water  molecules. 

FIG.  11.  Pair  energy  distribution  of  the  four  molecules  closest  to  the  HsOj  complex.  The 
dashed  line  is  for  the  classical  system  and  the  solid  line  is  for  the  quantum  system. 

FIG.  12.  The  rate  of  oxygen-oxygen  fluctuations  correlated  simultaneously  with  the  hydrogen 
bonding  of  a  water  molecule  adjacent  to  the  complex  (cf  Sec.  IV C).  The  quantum  rate  is  shown 
by  the  solid  line,  while  the  classical  rate  is  given  by  the  dash-dot  line. 

FIG.  13.  Power  spectrum  of  the  quantum  velocity  auto  correlation  function  of  the  O1H3  bond 
(solid  line).  The  broad  absorption  from  1000  cm“^  to  2500  cm“^  is  typical  of  strong  aqueous  acids. 
Shown  also  is  the  power  spectrum  of  one  of  the  covalent  OH  bonds  of  the  HsO^  complex  that 
does  not  participate  in  the  proton  transfer  (dashed  line). 
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